## $HeadUrl$
## $Id$
##=================================================
## Define NCounterSet object:

setClass('NCounterSet',
         contains='eSet');

## initializer
setMethod("initialize",
          signature(.Object="NCounterSet"),
          ## provide a signature consistent with the new defintion, but that picks up the (implicit) old
          function(.Object,
                   ##cdfName="", nrow=0, ncol=0,
                   phenoData = annotatedDataFrameFrom(assayData, byrow=FALSE),
                   featureData =  annotatedDataFrameFrom(assayData, byrow=TRUE),
                   experimentData=new("MIAME"),
                   annotation=character(0),
                   ##protocolData=phenoData[,integer(0)],
                   assayData=assayDataNew(
                     exprs = exprs,
                     ...),
                   exprs=matrix(numeric(0), nrow=nrow, ncol=ncol),
                   ## se.exprs
                   ...) {
            callNextMethod(.Object,
                           assayData=assayData,
                           phenoData=phenoData,
                           featureData=featureData,
                           experimentData=experimentData,
                           annotation=annotation)
          })

## accessor
setMethod("exprs",
          signature(object="NCounterSet"),
          function(object) assayData(object)[["exprs"]])

setReplaceMethod("exprs",
          signature(object="NCounterSet"),
          function(object, value) assayDataElementReplace(object, "exprs", value))

#if (is.null(getGeneric("exprs"))) setGeneric("exprs", function(object) standardGeneric("exprs"))
#if (is.null(getGeneric("exprs<-"))) setGeneric("exprs<-", function(object, value) standardGeneric("exprs<-"))

#setMethod("exprs", signature(object="NCounterSet"),
#          function(object) exprs(object))

#setReplaceMethod("exprs", signature(object="NCounterSet",value="matrix"),
#                 function(object, value) assayDataElementReplace(object, "betas", value))

panel.shade <-
function (x, y, ...) 
{
    r <- cor(x, y, use = "pair")
    ncol <- 14;
    pal <- col.corrgram(ncol)
    ##r.min <- min(r);
    ##r.max <- max(r);
    col.ind <- round(ncol * (r-.7)/.3);
    ##browser();
    usr <- par("usr")
    rect(usr[1], usr[3], usr[2], usr[4], col = pal[col.ind], 
        border = NA)
    rect(usr[1], usr[3], usr[2], usr[4], density = 5, angle = ifelse(r > 
        0, 45, 135), col = "white")
    box(col = "lightgray")
    text((usr[2]-usr[1])/2+usr[1], (usr[4]-usr[3])/2+usr[3], labels=round(r,2), col='white', cex=2)
}

panel.pts <-
function (x, y, pts.labels, ...) 
{
    plot.xy(xy.coords(x, y), type = "p", ...)
    box(col = "lightgray")
    text(x,y,labels=pts.labels, cex=.3, pos = 2);
}



### function that normalizes eset according to nCounter Data Analysis
### NCounterSet inherits from Biobase eSet

GetPosProbesExprs <- function(obj) {
  return(exprs(obj[which(fData(obj)[,"Code.Class"] == "Positive"),]));
}

GetNegProbesExprs <- function(obj) {
  return(exprs(obj[which(fData(obj)[,"Code.Class"] == "Negative"),]));
}

GetControlProbesExprs <- function(obj, control.probes) {
  ##browser();
  exprs <- exprs(obj);
  ##browser();
  return(exprs[match(control.probes, row.names(exprs)),]);
}

NormalizeNCounterNanostring <- function(obj) {
  ## normalize based on nanostring's suggestions

  ## average the positive probes
  exprs <- exprs(obj);
  pos.probes.exprs <- GetPosProbesExprs(obj);
  neg.probes.exprs <- GetNegProbesExprs(obj);
  
  pos.probes.exprs.colsums <- colSums(pos.probes.exprs);
  pos.probes.factor <- mean(pos.probes.exprs.colsums)/pos.probes.exprs.colsums;
  
  ret.exprs <- t(apply(exprs,1,function(x) {x*pos.probes.factor}));
  
  neg.probes.exprs.colmeans <- colMeans(neg.probes.exprs)*pos.probes.factor;
  ret.exprs <- ret.exprs - neg.probes.exprs.colmeans;
  
  ## set negative values to the average background values

  neg.indices <- which(ret.exprs <= 0);
  ret.exprs[neg.indices] <- neg.probes.exprs.colmeans[((neg.indices-1) %% ncol(neg.probes.exprs))+1];
  
  exprs(obj) <- ret.exprs;
  
  return(ret = obj);

}

NormalizeNCounterVSN <- function(obj, rep.samples = sampleNames(obj)) {
  ## normalize using VSN based on replicates.  by default, we assume that most lanes are about the same
  ## this is a DANGEROUS assumption since we only have ~100 genes

  ret.obj <- obj;
  
  norm.obj <- vsn2(exprs(obj[,rep.samples]), lts.quantile = 1);

  ## normalize to the replicates
  ret <- vsn2(exprs(obj), reference=norm.obj, lts.quantile=1);

  ##ret <- vsn2(exprs(obj));

  exprs(ret.obj) <- exprs(ret);
  
  return(ret.obj);
}

NormalizeNCounterVSN2 <- function(obj, control.probes) {
  ## normalize using VSN based on replicates.  by default, we assume that most lanes are about the same
  ## this is a DANGEROUS assumption since we only have ~100 genes

  ret.obj <- obj;
  
  norm.obj <- vsn2(exprs(obj[control.probes,]), lts.quantile = 1, minDataPointsPerStratum = 1, calib="none");

  ## normalize to the replicates
  ret <- vsn2(exprs(obj), reference=norm.obj, lts.quantile = 1, minDataPointsPerStratum = 1, calib="none");

  exprs(ret.obj) <- exprs(ret);
  
  return(ret.obj);
}

NormalizeNCounterQuantile <- function(obj) {
  ## normalize using quantile normalization
  ret.obj <- obj;
  
  norm.obj <- normalize.quantiles(exprs(obj));

  exprs(ret.obj) <- norm.obj;
  rownames(exprs(ret.obj)) <- rownames(exprs(obj));
  colnames(exprs(ret.obj)) <- colnames(exprs(obj));
  
  return(ret.obj);
}

### 
NormalizeNCounter <- function(obj, control.probes = featureNames(obj), rep.samples = sampleNames(obj), method = 'nanostring') {
  if (method == 'nanostring') {
    ret <- NormalizeNCounterNanostring(obj);

    if(!is.null(control.probes)) {
      ##print("here i am");
      control.probes.exprs <- GetControlProbesExprs(ret, control.probes);
      control.probes.exprs.colsums <- colSums(control.probes.exprs);
      
      control.probes.factor <- mean(control.probes.exprs.colsums)/control.probes.exprs.colsums;
      ##tapply(control.probes.exprs.colsums, as.factor(baseline.lanes), mean)[as.factor(baseline.lanes)]/control.probes.exprs.colsums;
      #print(control.probes.factor);
      ##browser();
      exprs(ret) <- t(apply(exprs(ret),1,function(x) {x*control.probes.factor}));
      ##colnames(ret.exprs) <- colnames(exprs);
      ##rownames(ret.exprs) <- rownames(exprs);
      ##browser();
    }    
  }

  if (method == 'log') {
    ret <- obj;
    exprs(obj)[which(exprs(obj)<=0)] = 1;
    exprs(ret) <- log(exprs(obj));
  }
  
  if (method == 'vsn') {
    ret <- NormalizeNCounterVSN(obj, rep.samples);
  }

  if (method == 'vsn2') {
    ret <- NormalizeNCounterVSN2(obj, control.probes);
  }

  if (method == 'quantile') {
    ret <- NormalizeNCounterQuantile(obj);
  }
  
    ##set negative values to 1
    ##ret.exprs[which(ret.exprs<0)] = 1;
    
    ##browser();
  return(ret);
}

CalculateFoldChange <- function(obj, baseline.lanes, log=TRUE) {
  #ret <- t(apply(exprs(obj),1,function(x) {x/exprs(obj)[,baseline.lanes]}));

  ret.obj <- obj;

  if(log == FALSE) {
    exprs(ret.obj) <- exprs(obj)-exprs(obj)[,baseline.lanes];
  } else {
    exprs(ret.obj) <- exprs(obj)*(1/exprs(obj)[,baseline.lanes]);
    exprs(ret.obj) <- log(exprs(ret.obj));
  }
  
  return(ret.obj);
}

CorrelateSamples <- function(obj) {
}

CombineDataSets <- function(objs) {
  ret.obj <- objs[[1]];
  
  for(obj in objs[2:length(objs)]) {
    ret.obj <- Biobase::combine(ret.obj, obj);
  }

  return(ret.obj);
}

PlotCorrelation <- function(obj, expr1, expr2, log=FALSE) {
  ##exprs(obj) <- exprs(obj)[-c(which(fData(obj)[,"Code.Class"] == "Positive"),
  ##                            which(fData(obj)[,"Code.Class"] == "Negative")),];
  
  if(log) {
    exprs <- log(exprs(obj));
  } else {
    exprs <- exprs(obj);
  }
  
  plot(exprs[,expr2],exprs[,expr1]);
  abline(lm(exprs[,expr1]~exprs[,expr2]-1));
}

PlotReplicates <- function(obj1, obj2) {
  plot(exprs(obj1), exprs(obj2), xlab = "Expression 1", ylab = "Expression 2");
  text(exprs(obj1), exprs(obj2), labels = featureNames(obj1), pos=1);
}

PlotVarianceVsMean <- function(obj) {
  plot(apply(exprs(obj),1,mean), apply(exprs(obj),1,sd), xlab = "Mean", ylab = "Variance");
}

PlotHeatmapExprs <- function(obj, log=FALSE) {
  ##exprs <- exprs(obj);
  if(log) {
    exprs(obj) <- log(exprs(obj));
  }

  heatmap.2(as.matrix(exprs(obj)), col=redgreen(70), trace = 'none');
}

PlotHeatmapCors <- function(obj, method = 'pearson', labels = NULL, pts.labels = NULL, ...) {
  ## first remove the positive and negative controls
  ##pos.probes.exprs <- GetPosProbesExprs(obj);
  ##neg.probes.exprs <- GetNegProbesExprs(obj);

  ##exprs(obj) <- exprs(obj)[-c(which(fData(obj)[,"Code.Class"] == "Positive"),
  ##                            which(fData(obj)[,"Code.Class"] == "Negative")),];

  ##if(log) {
  ##  exprs(obj) <- log(exprs(obj));
  ##}

  ## set negative infinite values to 0
  exprs(obj)[which(exprs(obj)<0)]=0;
  
  if(method == 'ccc') {
    cors <- CalculateCCC(obj);
  } else {
    cors <- cor(exprs(obj),use = 'pairwise.complete.obs', method=method);
  }
  ##cors[is.na(cors)] <- 1;

  corrgram(exprs(obj), order=NULL, labels=labels, panel=panel.shade, pts.labels=featureNames(obj), lower.panel=panel.shade, upper.panel=panel.pts, text.panel=panel.txt);
  
  ##heatmap.2(cors[cols,cols], cellnote=round(cors[cols,cols],digits=2), notecex=1, notecol='black',trace = 'none', cexRow=0.7, cexCol=0.7, na.rm=FALSE, dendrogram='none', margins = c(10,10), key=FALSE, Rowv=NULL, Colv=NULL, ...);
  ##heatmap.2(cors[cols,cols], trace = 'none', cexRow=0.7, cexCol=0.7, na.rm=FALSE, dendrogram='none', margins = c(15,15), key=FALSE, Rowv=NULL, Colv=NULL, ...);
  ##print(mean(cors));
}

CalculateCCC <- function(obj) {
  exprs <- exprs(obj);
  exprs.sd <- apply(exprs, 2, sd);
  exprs.mu <- apply(exprs, 2, mean);
  cors <- cor(exprs,use='pairwise.complete.obs', method='pearson');

  ccc <- matrix(data=NA, nrow=length(exprs.sd), ncol=length(exprs.sd), dimnames=list(colnames(exprs),colnames(exprs)));
  
  for(i in 1:length(exprs.sd)) {
    for(j in i:length(exprs.sd)) {
      ccc[i,j] <- 2*cors[i,j]*exprs.sd[i]*exprs.sd[j]/(exprs.sd[i]^2+exprs.sd[j]^2+(exprs.mu[i]-exprs.mu[j])^2);
      ###browser();
      ccc[j,i] <- ccc[i,j];
#      if(is.na(ccc[i,j])) {
#        browser();
#      }
    }
  }

#  browser();

  return(ccc);
}
